Two-dimensional Bloch electrons in perpendicular magnetic fields: an exact 
calculation of the Hofstadter butterfly spectrum 
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The problem of two-dimensional, independent electrons subject to a periodic potential and a 
uniform perpendicular magnetic field unveils surprisingly rich physics, as epitomized by the fractal 
energy spectrum known as Hofstadter's Butterfly. It has hitherto been addressed using various 
approximations rooted in either the strong potential or the strong field limiting cases. Here we 
report calculations of the full spectrum of the single-particle Schrodinger equation without further 
approximations. Our method is exact, up to numerical precision, for any combination of potential 
and uniform field strength. We first study a situation that corresponds to the strong potential limit, 
and compare the exact results to the predictions of a Hofstadter-like model. We then go on to 
analyze the evolution of the fractal spectrum from a Landau-like nearly-free electron system to the 
Hofstadter tight-binding limit by tuning the amplitude of the modulation potential. 



The motion of electrons in a crystalline solid subject 
to a magnetic field has been considered since the early 
days of quantum mechanics [1 j. The field splits the crys- 
tal's electronic bands into sub-bands and internal mini- 
gaps; the energetic arrangement of these sub-bands forms 
a fractal structure reminiscent of a butterfly when plot- 
ted as a function of the field [8]. Significant experimen- 
tal effort has been devoted to detecting signatures of 
this energy spectrum in two-dimensional electron gases 
(2DEGs) [2j|3]. For independent electrons, the situation 
is described by the single-particle Schrodinger equation, 



Hi/j(r) 



1 

2m J 



V(r) 



^(r) = £ty(r), (1) 



where H is the Hamiltonian, ^(r) is an eigenstate with 
energy E, II = p + eA(r) is the dynamical momentum 
operator and A(r) is the vector potential corresponding 
to the magnetic field, B = V x A. We take the field to 
be uniform and oriented along the ^-direction, B = B e z . 
The electrons are restricted to the two-dimensional (2D) 
xy-plane, and the external potential V(r) is periodic on 
a Bravais lattice defined by vectors 



R n =ja + £;b, n = (j,fc)GZ 2 . 
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To date, this problem has been chiefly approached by 
approximations starting from two complementary limits, 
considering either the influence of a weak magnetic field 
on the band structure resulting from a strongly varying 
potential ^(r) [HE], or the influence of a small mod- 
ulation potential on the Landau-quantized electrons in 
a strong field [5j [6]. In the strong potential limit, one 
typically starts with a tight-binding (TB) approxima- 
tion for a single band of the zero-field (A = 0) prob- 
lem, -E'(k), where k is the crystal momentum. Then, 



an effective Hamiltonian for the magnetic field problem 
is generated from E(k) through the Peierls substitution 
k — > (p + eA) /h [TJE]. This procedure was used, among 
others, by Hofstadter in his seminal article for a nearest- 
neighbor (NN) TB model of the 2D square lattice [8]. 
There are a number of simplifications inherent to this ap- 
proach: (i) the TB approximation of the zero-field band 
structure, (ii) the restriction to electrons in a single band, 
(iii) the neglect of the diamagnetic energy of the TB or- 
bitals and the field dependence of the TB hopping in- 
tegrals. This has been shown to lead to quantitative as 
well as qualitative errors for both nearly free and tightly 
bound two-dimensional electrons [5HT2]. Generalization 
of the effective Hamiltonian approach has turned out to 
be difficult, see, e.g., Ref. [13 . Surprisingly, the strong 
field approach is closely related to the strong potential 
one: if potential-induced coupling between different Lan- 
dau levels is neglected, the same secular equation results, 
but with the magnetic field replaced by its inverse [5]. 
Including such coupling has a profound effect on the cal- 
culated energy spectrum p~4j [15] , even for weak coupling 
strength. The resulting rearrangement of the Hofstadter 
butterfly has recently been confirmed by experiments [3] ; 
we take this as a strong indication that experiments can 
only be fully understood by going beyond the approxi- 
mate schemes described above. 

Method. In the absence of magnetic field, Bloch's the- 
orem, which results from the commutation of the lattice 
translations with the Hamiltonian, allows to restrict the 
calculation of the eigenfunctions ^(r) to one primitive 
cell of the lattice in Eq. As these translations do not 
leave the vector potential A(r) invariant, they no longer 
commute with H when a magnetic field is present; con- 
sequently, the eigenfunctions ^(r) are not Bloch waves. 
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By combining a lattice translation with a suitable gauge 
transformation to counteract its effect upon the vector 
potential, one can define magnetic translation operators, 
^a(Rti)? that do commute with H, see Refs. [T6HT8] . The 
operators T^(R n ) do not, however, form a group. For a 
rational field, where the number of magnetic flux quanta 
per unit cell, a = ^| (ax b)5, is a rational number, 
a = p/q with p and q relative prime, one can choose a 
larger magnetic lattice that has an integer number of p 
flux quanta passing through each cell, e.g., 



S n =ja + %b), n=(j,k)eZ 2 
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On this lattice, a generalized version of the Bloch theo- 
rem holds [18]: for an orthorhombic lattice and Landau 
gauge, A(r) = Bxe y , the solutions of Eq. can be 
chosen of the form </>(r) = e l6r u e (r), where u e (r) obeys 



u e (r + S) =exp 
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and the magnetic crystal momentum is restricted to 
the first Brillouin zone of pi). This condition allows to 
restrict the calculation of %r (r) to one primitive cell of 
the magnetic lattice. The functions u e (r) exhibit a pe- 
culiar topology in a space where their ^-component is 
expanded in a Fourier series, see Ref. [19] . This allows to 
take into account the boundary condition Q in a natu- 
ral way when the diffusion method is used for solving the 
Schrodinger equation [20, 21]. We have summarized the 
steps leading to Eq. Q and the technical details of the 
diffusion method in the Supplementary Information. 

At a fixed rational field, a = p/q, we obtain the bands 
Ej (0) by numerically solving Eq. ([I]) for a grid of 0- values 
spanning the magnetic Brillouin zone (MBZ) correspond- 
ing to this field. The density of states (DOS) as a function 
of the field, p(B,E), can then be calculated by integrat- 
ing the bands over the MBZ and repeating the process 
for different fields. The rational field can only be tuned 
in discrete steps, ap = P/Q, where P, Q G Z. The size of 
the magnetic unit cell depends on the reduction of P/Q 
to a quotient of relatively prime integers p/q, and is q 
times as large as the zero-field ("geometric") unit cell. 
At the field ap, one zero-field band is thus expected to 
split into q magnetic bands, which are known to cluster 
in groups of p bands [8]. The pattern of distinct prime 
factors of P and Q as ap is swept across a range of fields 
gives rise to the self-similar, fractal pattern of gaps in the 
DOS that has become known as the Hofstadter butterfly. 

When the Fermi energy of the system lies in a gap, i.e., 
a region where p(B, Ep) is zero, the Hall conductance a xy 
assumes a quantized value, <rf^ p = ne 2 /h,n G Z. The 
fundamental topological reason for this quantization was 
revealed by Thouless et al [22], who showed that both 
in the strong field and strong potential limits the Kubo- 
Greenwood formula for <j xy is related to the Chern num- 
ber of the U(l) bundle over the magnetic Brillouin zone. 
This was later argued [23 to be a direct consequence of 



the magnetic translation symmetry, Eq. Q • Sweeping ei- 
ther the magnetic field or the Fermi energy through the 
fractal pattern of minigaps inside a broadened band or 
Landau level results in a peculiar, no n- monotonous Hall 
effect [22 . Indications of this behavior have been found 
in experiments [2 . In this work we have used an alterna- 
tive approach, introduced by Stfeda [24], to obtain crf^ p 
from the numerically calculated DOS, 



a^(B,E F ) = e 
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In our numerical calculations, we have studied a sim- 
ple system consisting of a two-dimensional square lattice 
of potential wells with the symmetrized Fermi function 
form [25] , 



^) = t/cothg 



cosh (§) + cosh (a) 



with parameters 



r = 39.7 nm, d = 1.59 nm, a = 100 nm, 



(6) 



(7) 



where a is the lattice spacing. The potential is illustrated 
in Fig. ([l^); the parameters are set to loosely reproduce 
the conditions of earlier experimental studies [2j [3] . 

Strong Potential Regime. We first compare the spec- 
trum generated by the full Schrodinger equation ([!]) to 
the results of a TB approximation similar to the one used 
by Hofstadter. We thus choose a fairly deep modulation 
potential, 



U = -V = -8.4 meV. 



(8) 



We obtained the band structure and DOS at zero mag- 
netic field [see Fig. ([IJ)] by numerically solving the corre- 
sponding single-particle Schrodinger equation. We then 
fitted a TB model to the lowest three bands [black dashed 
lines in Fig. ([IJ)], and employed the Peierls substitution 
to obtain the fractal energy spectra shown in Fig. ([l]i), 
see the Supplementary Information for details. The low- 
est band [n = 1 in Fig. ([1J)] required only nearest- 
neighbor hopping integrals in the TB band model for an 
adequate fit, and thus yields a spectrum corresponding to 
that obtained by Hofstadter [8] [left panel of Fig. ([l]i)]. 
The second and third bands required up to 3 rd -nearest- 
neighbor hoppings, which lead to significantly distorted 
versions of Hofstadter's butterfly, shown for the second 
band on the right panel of Fig. ([l]i). The spectrum of 
the third band is similar and is not shown. Our findings 
qualitatively agree with the results of Ref. [26] . 

In Fig. ([2^l) we show the energy spectrum of the six 
lowest bands, as obtained by numerically solving the full 
magnetic eigenvalue problem, Eq. 0, using the scheme 
outlined above. A maximum magnetic unit cell size of 
Q = 32 was employed in the calculation. For compari- 
son, the spectrum of an isolated potential well is plotted 
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FIG. 1: (color online), a) Schematic plot of the periodic 
Fermi well potential, b) First Brillouin zone of the recip- 
rocal lattice, c) Left panel: Zero-field band structure of the 
potential with parameters |7| and U — — Vq. Solid red lines 
show bands obtained by numerically solving the Schrodinger 
equation, dashed black lines are TB bands fitted to these ex- 
act bands (see text). Right panel: DOS p(E) obtained from 
the full Schrodinger equation, d) Hofstadter butterflies for 
the two lowest bands. Areas with non-zero DOS are printed 
black, and the gaps are colored according to the correspond- 
ing quantized Hall conductance crf^ p in units of e 2 /h. White 
indicates zero Hall conductance, warm (cold) colors indicate 
positive (negative) Hall conductance. A number of larger gaps 
are labeled with the corresponding Hall conductance for ref- 
erence. The butterflies are periodic in the flux, one period 
being shown in each case. 



in Fig. (|2p). It bears a strong resemblance to the Fock- 
Darwin (FD) spectrum [27] of a parabolic well, we will 
thus refer to these states as "FD states" in the following. 
It can be seen in Fig. that in the periodic system 
the FD states of the isolated well are broadened into 
bands with a fractal internal structure that is qualita- 
tively well described by the Hofstadter butterfly. In the 
exact result the periodicity of the Hofstadter spectrum 
is superimposed onto the field dependence of the corre- 
sponding FD state; this field-dependence is not taken into 
account in the TB model with constant hopping integrals. 
In general, higher energy levels, having more extended 
wave functions, undergo larger broadening at a given flux 
value. Conversely, bands become narrower with increas- 
ing field, as their wave functions become more spatially 
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FIG. 2: (color online). Magnetic energy spectrum (DOS and 
Hall conductance <r|^ p ) for the lowest six bands of the of the 
square Fermi well lattice with parameters Q and U = — Vb, 
calculated by numerically solving Eq. |l]). The color coding 
and labeling is as in Fig. |l]). Panel (b) shows a magnified por- 
tion of the lowest band indicated by the dashed box in panel 
(a). Panel (c) shows the energy spectrum of an isolated Fermi 
well potential, Eq. §6§. It strongly resembles the spectrum of 
a parabolic potential (the Fock-Darwin spectrum [27]). 



localized, tending to Landau levels in the limit of high 
field intensities. In the TB model, the main effect of in- 
creasing the second- and third-nearest neighbor hoppings 
is a distortion of the butterfly that opens a gap at flux 
strengths of ctj = j + 1/2, with j £ Z, indicated by the 
vertical (red) dotted lines in the butterflies in Fig. ([l]i). 
This behavior is also present in the exact results, compare 
regions (4)-(5) or (7)-(8) in Fig. The gap widens 

for higher-energy bands, consistent with the TB approx- 
imation, where such bands need to be modeled by larger 
hopping integrals to more distant neighbors. The gap 
decreases again with increasing field due to the stronger 
localization of the wave functions, an effect that is not 
included in the TB description with ^-independent hop- 
ping. Band crossings, which are not described within 
the single-band Peierls approximation, are an interesting 
subject for future studies: the narrow third FD band 
seems to disrupt the butterfly pattern of the broader 
bands it crosses in regions (3) and (6), but does not seem 
to exert any noticeable effect in region (9) . At higher en- 
ergies [see region (10)], where multiple bands cross, the 
resulting fractal spectrum can assume a form that is very 
different from the original Hofstadter butterfly. 
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Intermediate Potential Regime. We now explore the 
evolution of the spectrum as the potential is changed 
from nearly flat (strong field limit) to highly modulated 
(strong potential limit). The intermediate stages of this 
evolution are not accessible to the approximate method- 
ologies hitherto employed. We use the square lattice po- 
tential with the parameters given in Eq. (J7|, but change 
the well depth U. We start with a very shallow potential 
(U = — 0.06 Uo), which corresponds to a Landau-like sys- 
tem with nearly-free electrons in a magnetic field. The 
resulting spectrum is illustrated in Fig. (|3^i); it exhibits 
the typical "Landau fan" form, with slightly broadened 
Landau levels (LLs) that display an internal fractal struc- 
ture of minigaps (most evident in the lowest level, Lo). 
The Hall conductance between LLs increases monotoni- 
cally in steps of e 2 //i, consistently with the integer quan- 
tum Hall effect. Furthermore, an emerging white gap 
( a xy P ~ 0) can be observed at low field (a = — 1, 
E « 0.11 Vb). When the potential modulation is in- 
creased to U = — 0.13Vb [Fig. ^jp)} this gap is seen to 
widen further, practically cutting off a low-energy trian- 
gular section (marked "T" in the figure) from all LLs 
L n with n > 0. At a = 1 the upper-right tip of this 
triangular cluster of states retains a tenuous link to the 
broadened L\ Landau band, at the position indicated by 
the blue arrow. At the same time, the minigaps at a > 1 
in the lowest Landau band Lo have broadened further, to 
the point where the band is only held together by a nar- 
row sub-band at a = 1, indicated by the red arrow. The 
Hall conductance in the gap between Lo and L\ at a = 1 
is still <rf^ p = 1. But upon further increasing the modu- 
lation strength, the gap first closes and then reopens with 
the above links reversed: at U = — 0.26Uo [Fig. ([3^)], the 
tip of the triangular cluster is now connected to the low- 
est miniband of Lo (blue arrow), while the rest of the Lo 
miniband is connected to L\ (red arrow). At the same 
time, the gap has changed its character to <rf^ p = 0, re- 
sulting in the first clearly discernible FD-like band Fo, 
separated from all higher energy states by an unbroken 
white gap. Most of Fo is formed from the lowest mini- 
band of the LL Lo, except for the triangular cluster of 
states in a < 1, which originates from higher LLs. The 
process of formation of the first FD band (Fo) that is 
incipient at a = 1 in Fig. ^$jp) can be seen to repeat 
itself at a = 2 for the second FD band (F x ) [Fig. 
green and purple arrows] and at a = 3 for the third FD 
band (F2) [Fig. (|3ji)] - The gradual transformation from 
LLs to FD-like bands thus proceeds by the re-connection 
of minibands from one broadened LL to those of neigh- 
boring LLs. This rearrangement process fragments the 
large triangular gaps between LLs in the weak potential 
limit into minigaps encapsulated by the emerging FD-like 
bands. As a result, the FD states F , Fi, F 2 in Fig. [3]i 
are composed of a low-field section that originates from 
higher LLs (triangle T in case of F , one Hofstadter but- 
terfly segment plus a triangular cluster in case of Fi, 



etc.), and a high-field section that is one miniband of the 
lowest LL Fq. 

The above results show that with the method pre- 
sented here it is possible to perform exact calculations 
of the spectrum of independent electrons in a 2D pe- 
riodic potential and constant perpendicular magnetic 
field. The same technique is readily applicable to solve 
the Kohn-Sham equations of density functional the- 
ory (DFT) which are expected to provide a reasonable de- 
scription at least for weakly correlated electron systems. 
Experimental techniques have recently become available 
to directly probe the local DOS in 2DEGs on surfaces 
in a perpendicular magnetic field using scanning tunnel- 
ing spectroscopy, making possible the detection of spa- 
tial features of some LLs [32 , and even measuring their 
response to a ID periodic potential due to surface buck- 
ling [33]. Such developments may presently bring about 
additional possibilities of experimentally determining the 
spectral features of the 2DEG subject simultaneously to 
a periodic potential and a perpendicular magnetic field, 
contrasting them with the predictions reported herein. 
The Hofstadter butterfly is not unique to the system dis- 
cussed here. The basic principle underlying the fractal 
energy spectrum is the presence of two competing sym- 
metries (in the case studied here, the periodicity of the 
lattice and the symmetry of the Landau orbits, which is 
governed by the area required by one magnetic flux quan- 
tum). Similar patterns have been observed or predicted 
to occur in a variety of very different systems, such as 
microwaves transmitted through a waveguide with a pe- 
riodic arrangement of scatterers [28] . the electronic [29] 
and vibrational [30] spectra of incommensurate crystals, 
ultracold atoms in optical lattices [31] , and photonic crys- 
tals [34] [35] . The topological protection of the quantum 
Hall phase has been shown to improve the performance 
of optical delay lines and to overcome limitations related 
to disorder in photonic technologies [35] . We hope that 
the method presented here will also prove beneficial in 
these related fields. 
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